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ABSTRACT: We show that the dynamics of a spatially closed Priedmann - Robertson 
- Walker Universe conformally coupled to a real, free, massive scalar field, is chaotic, for 
large enough field amplitudes. We do so by proving that this system is integrable under 
the adiabatic approximation, but that the corresponding KAM tori break up when non 
adiabatic terms are considered. This finding is confirmed by numerical evaluation of the 
Lyapunov exponents associated with the system, among other criteria. Chaos sets strong 
limitations to our ability to predict the value of the field at the Big Crunch, from its given 
value at the Big Bang. 
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I - Introduction 

In this paper, wc shall present analytical and numerical evidence of chaotic behavior 
in a cosmological model, consisting of a spatially closed, Friedmann - Robertson - Walker 
(FRW) Universe, filled with a conformally coupled but massive real scalar field, although 
this model is far too simplified to be considered realistic, its simplicity itself makes it 
an interesting testing ground for the implications of chaos for cosmology, either classical, 
semi - classical or quantum. Moreover, the occurrence of chaotic motion in this example 
supports the conjecture of a farther reaching connection between chaos at the classical level, 
particle creation at the semiclassical one, and decoherence in the full quantum cosmological 
treatment ( see below ). 

The research on chaotic cosmological models began with the work by the Russian 
school'^] and by C. Misner and coworkers on chaos in Bianchi type IX vacuum cos- 
mologies, although the hopes initially placed on these models ( such as a resolution of 
the horizon problem from the mixing effect associated to chaos ) could not be sustained, 
Bianchi IX is still by far the best studied chaotic cosmology , both because of its intrin- 
sic interest, and as a test bench for more general questions, such as whether chaos is an 
observer dependent phenomenon t^l . In the broader context of General Relativity, chaos 
has been considered in connection to geodesic motion, both in cosmological and Black 
Hole space times t^l. 

The relative paucity of examples of relativistic chaos^^l makes it hard to disentangle 
the general features ( if any ) of these phenomena from the "miracles" proper to each 
peculiar manifestation, such as the occurrence of the Gauss' map hidden in Bianchi IX 
dynamics, which led to the discovery of cosmological chaos in the first placet-*^! . It is there- 
fore of the utmost importance to develop a systematic search for instances of relativistic 
chaos, to enhance our battery of examples, and therefore better to aim future research. 

As pointed out in an earlier communication , the class of near integrable systems 
(NIS) is an interesting field for such a search. A NIS is an integrable system which becomes 
non integrable under the effect of a perturbation. The non integrable perturbed dynamics 



2 



may or may not be, in turn, chaotic. A typical situation occurs when the unperturbed 
dynamics displays an unstable fixed point, asymptotically joined to itself by a nontrivial 
orbit ( the so called "homoclinic loop" ). If the loop is destroyed by the perturbation, then 
a "stochastic layer" forms in its neighborhood. There are subsets in this layer where the 
dynamics is equivalent to a Bernoulli shift . This is the same degree of chaotic behavior 
characteristic of the Bianchi IX example. A detailed discusion of this scenario, geared to 
its application in relativistic problems, has been given in Ref.8. 

The case at hand in this paper belongs to a wider class where there are no homoclinic 
points in the unperturbed dynamics. However, in a certain sense, the perturbation both 
creates them and destroys the corresponding homoclinic loops, allowing the formation of 
stochastic layers. This road to chaos results from the overlapping of several resonances 
between the external perturbation and the unperturbed motion. In the neighborhood of a 
resonant region in phase space, one of these resonances dominates over the rest; because 
of this dominant resonance, KAM tori are destroyed in this region. In particular, the 
torus for which the resonance is exact dissapears, leaving behind a discrete set of both 
stable and unstable fixed points'-'^^'. The unstable fixed points are connected to each other 
through doubly asymptotic orbits or separatrices; chaos occurs whenever the separatrices 
are further destroyed by the effect of the secondary resonant terms. In the neighborhood 
of a destroyed separatrix, it is possible to find ( measure zero ) invariant sets where the 
restricted dynamics is equivalent to a Bernoulli flow. More generally, chaotic behavior, 
under the form of extreme sensitivity regarding initial conditions, is displayed in a set of 
full measure around the separatrix, the so called stochastic layer t^^L 

The particular model we shall present here has been chosen both because of its 
simplicity and because of its relevance to the discussion of the ties between classical and 
quantum cosmological modelst^^L Indeed, similar models have been used by Hawking and 
Paget^"^] to discuss the relationship between the cosmological and thermodynamic arrow of 
time, in the framework of Quantum Cosmology. Since in this model the system behaves as 
a Bernoulli flow in appropiate invariant sets, time reversal invariance may be spontaneously 
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broken, and thus a "chaotic" arrow of time appearst^^^. The behavior of this arrow supports 
Page's picture of the relationship between the cosmological and the thermodynamic ones. 

In the cosmological literature more generally, fundamental scalar fields are usually 
considered within the framework of Inflationary models. In these models, however, minimal 
coupling is often preferred to conformal onet^^L Thorough analysis of the dynamics of a 
minimally coupled field in a FRW background have been developped by several authorst^^l; 
they have found that Infiationary cosmologies act as an attractor in phase space, but little 
or no evidence of chaos. Further studies of non minimally coupled fields in open universes 
have yielded similar resultst-*^^!. In the conformally coupled model we shall present here, 
on the other hand, "inflationary" periods, where the radius a of the Universe increases by 
several orders of magnitude, are a common feature of the solutions to the full dynamics. 
In this model, "inflation" is not powered by an effective cosmological constant, but rather 
by the effective Newton's constant becoming negative for large values of the field '^^1 . 

Another source of interest in this model is that it is one of the simplest cosmologies 
where particle creation by the gravitational field occurs t^^l. Indeed, particle creation in 
models of this kind has been analyzed by Hartle and others '^^l . In a perturbative analysis, 
the spectrum of created particles is closely related to the Fourier decomposition of 0^(77), 
where r] stands for conformal time, and a is the "Radius of the Universe" ( see below ). 
This fact, and the connection between particle creation and the breakdown of the WKB 
approximation through Stokes' phenomenon t^^l, suggest that particle creation is indeed 
the effect of resonances between the evolution of a and that of the scalar field. If this were 
true, then a strong correlation between classical chaos and semiclassical particle creation 
should exist. The relationship of semiclassical particle creation to quantum cosmological 
decoherence has been discussed elsewhere t^^^. 

In this paper we shall investigate our cosmological model both through analytical 
(perturbative) and numerical (non perturbative) methods. After showing how the full 
Hamiltonian can be split into an integrable part plus a perturbation, we shall show how 
the KAM tori of the integrable part are destroyed and replaced by stochastic layers under 
the effect of the perturbation. We shall analyze the onset of chaos through three methods 
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of increasing sophistication. We shall begin by employing Chirikov's resonance overlap 
criteriumt-'^^l to show that the stochastic layer of a given resonance is wide enough to in- 
trude into the layers of its neighbors. Then we shall employ the Melnikov Integral^^^ to 
show that the separatrix created by each resonance is split by the effect of the secondary 
perturbations. Finally, we shall investigate the dynamics close to a separatrix by construct- 
ing the corresponding separatrix or standard map ^^^^ We shall use this map to discuss 
the important issue of whether the sensitive dependence on initial conditions proper to 
chaos is strong enough to lead to observable effects within the lifespan of the Universe. It 
should be remembered that Bianchi IX models have been found lacking in this regard^^^ . 

Numerical analysis will allow us to go beyond perturbation theory. We shall present 
results from numerical integration of the full equations of motion, both in the chaotic 
and non chaotic regimes; plots of the field against a, where the change in the topology 
of the orbits, subsequent to chaos, can be observed; numerical estimates of the Lyapunov 
coefficients [^^1 in the stochastic layers; and finally, a plot of the values of the field at the 
Big Crunch vs. those at the Big Bang, for fixed initial geometry, to demonstrate the loss 
of predictibility subsequent to chaos. 

The paper is organized as follows. We introduce the model in the following section, 
where we also develop its perturbative analysis. Our numerical results are presented in 
Section III. Finally, we briefly state our conclusions in Section IV. 

II - The model and its perturbative treatment 

III - The model 

Let us begin by introducing the model and how we shall split its Hamiltonian into 
an unperturbed part plus a perturbation. 

As already discussed, our cosmological model assumes a FRW spatially closed geom- 
etry, that is, a line element 
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(1) 



Where < 99 < 27r, O<0<7r, 0<x<7r, and rj stands for "conformal" time. 
For concreteness, we shall consider only models starting from a cosmic singularity, that is, 
we restrict r] to be positive, with a(0) = 0. Also, as indicated by the dynamics, we shall 
assume that after the Big Crunch ( that is, when a returns to ), a new cosmological 
cycle begins, now with a < 0. Therefore, a complete periodic orbit describes the birth and 
death of two Universes. 

The gravitational dynamics is described by the Einstein - Hilbert action 




Where the determinant of the metric —g = a^sin^xsin^, and the scalar curvature 




(We shall use MTW conventions throughout t^^^) A dot represents a rj derivative, 
and mp is Planck's mass. For simplicity, we shall assume rup = •\/l/12v, where v = 27r^ is 
the conformal volume of an spatial surface. 

The action for a conformally coupled, massive, real scalar field is given by 

Sf^ J d'^x + (m^ + (l/6)i?)$2] (4) 

Where m? is the mass. Consistency with the symmetries of the background geometry 
demands the field be homogeneous. Parametrizing the field as $ = (/)/t'*^^/^-*a, performing 
the spatial integrals, and discarding total derivatives, we are led to a dynamical system 
with two degrees of freedom, a and 0, and Hamiltonian 
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H = ( J)[-(7r2 + a') + ip' + 4>') + m^a^cl>^] 



(5) 



Where tt and p are the momenta conjugated to a and respectively. At this point 
we must recall that, because we have relinquished our gauge freedom in writing the line 
element as in Eq.(l), we are missing one of Einstein's equations, namely, the Hamiltonian 
constraint [^^1 . We reintroduce this constraint as a restriction on allowable initial conditions 



Which is clearly respected by the dynamics. 

When VP? = 0, the Hamiltonian (5) is obviously integrable; for nonvanishing m^, this 
is no longer so obvious, and we must resort to perturbative methods ( or solve the dynamics 
numerically, see next section ). It is tempting to consider the massless Hamiltonian as 
the unperturbed one, with the last term in Eq. (5) as perturbation. This is questionable, 
however, on the grounds that for "macroscopic" Universes, very easily we obtain vr?o? ^ 1. 
For example, for our own Universe a ~ 10^°, and for a mass of 1 eV, we get m ~ 10~^®, 
so ma ~ 10^^. In this regime, the last term in (5) is in no way small compared with the 
other ones. 

We shall therefore proceed in a different way. Let us observe that, if the evolution of 
a is slow compared to the oscillations in then its main effect is to produce an adiabatic 
change in the frequency of the latter. We therefore introduce the "adiabatic" amplitude 
and phase j and <^ 



= 



(6) 




(7) 



p = -s/^ujj cos (f 



(8) 
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Where ijp" = l+m^a^ is the instantaneous frequency of the field. This transformation 
is canonical; it can be accomplished by means of the generating functional 



Si = Pa+{-^) (9) 

^2tan^^ ^ ^ 

But, because of the a dependence in a;, we are forced to change the geometrical 
momentum from tt to P, according to 



7r = P+——-^sm2cp 10 
2(1 + m^a^) 

We now substitute the new variables in the Hamiltonian, and rewrite this as 



H = -{Ho + SH) (11) 
where the unperturbed Hamiltonian 

^0 = {l)[P^ + a'] - jVl + m^a^ (12) 

is obviously integrable ( Hq and j are constants of motion in involution ), and the 
perturbation 



5H=——^ f^sin2(^+f ,^ , l-cos4yp 13 

2{l-\- m?a^) ^4(1 + m2a2)^ \ -r; \ / 

remains small both for small and large Universes. 



11.2 - Solving the unperturbed dynamics 

The dynamics of a, as generated by Hq, is obviously bounded. The point o = is a 
fixed point; it is stable if m^j < 1, , and unstable otherwise. In this second case, there is an 
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homoclinic loop associated with it. However, this orbit does not satisfy the Hamiltonian 
constraint, Eq. (6); rather, we have Hq = —j on the homochnic loop. 

The equations of motion are simpler if written in terms of a new variable X = 
Vl + rmP'o?-, rather than a itself. The transformation is canonical if we associate to X the 
momentum 



Px = (14) 
In terms of the new variables {X, Px), the unperturbed Hamiltonian reads 



Ho = r i P'x + - m'jf - {m'jf - 1] (15) 

For fixed Hq = h and j, provided h > —j, X ranges from 1 to the classical turning 
point Xt = m^j+my/K, where K = (l/m?) + (mj)'^ + 2h ( The reasons for the somewhat 
unusual notation shall be clear below). 

Although the Hamiltonian Eq. (15) can be explicitly integrated ( the solution involves 
the use of elliptic integrals '^^1), for our purposes it will be better to concentrate on a 
particular case, where a number of simplifications will be available. Concretely, we shall 
consider the large j limit, with Hq ~ and m fixed. In this limit, we find X ^ 1 for 
most of the orbit. Indeed, from Hamilton's equations, we get X ~ 2my/h + jy/X — 1 when 
X ^ 1 and Hq = h. Therefore, even if X starts at 1, X — 1 becomes of order unity after a 
lapse 5t ~ {my/h + which is small in the case under consideration. Now, for X ^ 1, 
the Hamiltonian Eq. (15) simplifies to 

Ho = (^)Pl + {^)[{X - m^jf - {m'jf - 1] (16) 
If we parametrize 



Px 



cos2q! 



m 
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(17) 



X = rn^j + my/K sin 2a 



(18) 



Eq. (16) reduces to 



Ho = {l)[K - imjf - ^] (19) 

As X goes from 1 to Xt, that is, over a quarter of an orbit, sin 2a goes from —{m^j — 
l)/my/K to 1. For large j and i?o ~ 0, however, K is very close to (mj)'^. So the end 
points can be taken as o; ~ ±7r/4, in which case K is precisely the associated action 
variable. We have succeeded in integrating the unperturbed motion, only that, because 
the reparametrization eqs. (17) and (18) involves j, the angle canonically conjugated to it 
is no longer (f, but 



e = ip- m\fK cos 2a (20) 

The transformation eqs. (17), (18) and (20), from variables {X, Px,(p,j) to new 
variables {a,K,9,j' = j) is canonical; indeed, it is generated by 



In the following, we shall omit the prime on j' . 



II. 3 - Analysis of the perturbation 

Having reduced the unperturbed Hamiltonian Eq. (12) to action - angle form (19), 
let us return to the perturbation, Eq. (13). We notice, first of all, that the first term 
in the perturbation dominates the second ( except close to the turning points ). For our 
purposes, it shall be sufficient to retain only the dominant perturbation. 

We also notice that, my/K being large, sin2(fi = sin2(^ + my/K cos2a) (cfr. Eq. 
(20)) is a strongly oscillating function of a. Thus, in computing the Fourier expansion of 
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the perturbation, we may substitute the non oscillatory factor m^aPj /2{1 + m^a^) by its 
mean value over a quarter orbit, which is easily found to be (j/tt) [InX^ — {X\ — l)/2Xy] ~ 
(jlnXT)/7r. 

The Fourier expansion of the oscillatory factor itself is 

oo 

sin 2{Q + mVK cos 2a) = Jo(2m\/K) sin 2^ + ^ Jn{2my/K)[smS+ + sin 5"] (22) 

n=l 

Where Jn is the usual Bessel function, and 5^ = {29 ± 2n{a + (n/A))). Therefore, 
resonances occur whenever 

Uj ± nujK = (23) 

^j,K being the frequencies associated to 9 and a, respectively. Prom Eq.(19), we have 
u>j ~ rri^j ( we have reinstated the proper sign, cfr Eq.(ll) ), and ujk ~ ~l/2- We thus find 
a tower of resonances, corresponding to all positive values of n; the n-th resonance occurs 
cit jn ~ n/2m^, independently of K. If we further impose the Hamiltonian constraint, 
then K must take the value Kn ~ (rnjn)^ = (n/2m)^. 

The n-th resonant term in the Hamiltonian, with its proper sign, reads 

( — )(lnXT) Jn(2myK) sin(2^ + 2n{a + (7r/4))) (24) 

TT 

To analyze the perturbed motion it is sufficient to approximate the prefactor by 
its value at resonance. So doing, and using the proper asymptotic form for the Bessel 
function'^^1 , Eq.(24) reduces to 

(^)n2/3(lnn) sin(26' + 2n{a + (7r/4))) (25) 
where e is a numerical coefficient, e ~ 0.111827... 
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11-4 - Solving the perturbed motion 

To analyze the motion in the presence of the perturbation, let us focus first on the 
neighborhood of the n-th resonance ( the calculations below are simplest if n = 4fc + 1 
for some k, which we shall assume ). Let us keep only the dominant resonant term Eq. 
(25) in the Hamiltonian, and introduce the new action variables ^ = (j — jVi)/2 and 
K, = —{K — Kn — n{j — j„))/2, which vanish at the resonant point. These variables are 
canonically conjugated to the angles ip = 2(0 + na) and r = —2a, respectively ( the 
transformation is generated by ^3 = [Kn — 2{k. — n$,)]a + [jn + 2$]9 ) . In terms of the new 
variables, the resonant Hamiltonian becomes ( cfr. Eqs. (19) and (25) ) 

iy„ = « + 2771^^2 -(^) cos V' (26) 
Where loq = 4en^/^ In n. 

Because the Hamiltonian (26) is linear in k, it can be considered as resulting from the 
"parametrization" of a one degree of freedom system'^^^ . To analyze the resulting motion, 
it is convenient to "deparametrize" it, that is, to promote r to the roll of "time". The 
evolution of the "true" degree of freedom '0 as r unfolds is described by the Hamiltonian 

E = -At. 

From Eq. (26) and the Hamiltonian constraint (6), it is obvious that E is simply the 
Hamiltonian of a non linear pendulum, being the frequency of small oscillations around 
the stable equilibrium point '0 = 0. The pendulum also has an unstable equilibrium point 
i/j = ±7r, joined to itself by a separatrix. 

Following Chirikovl^^l, we define the "width" of the resonance as the maximum value 
of ^ along the separatrix, ~ uJo/2m'^. Since uq ~ n^/^, it is clear that for large n each 
resonance is much wider that the separation l/2m? between resonances. Thus, according 
to Chirikov's criterium, the behavior of the perturbed system must be chaotic. 

To obtain a more detailed picture, let us add to the resonant Hamiltonian Eq. (26) 
the first secondary resonance (n = 4/c + 2), 
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6Hn ~ {coq/ 4:771^) sm('0 - r) 



(27) 



( The same term is added to the "deparametrized" Hamiltonian E). We can now 
apply the criterium that chaos will occur if the secondary resonance is able to destroy the 
separatrix created by the primary resonance. This can be determined by computing the 
total change in the unperturbed pendulum Hamiltonian, induced by the perturbation, as 
the pendulum swings along the separatrix, which is given by the Melnikov Integral'^! 

/+00 
(ir ^cos(V'-t) (28) 
-oo 

In our case'-'^-'^^, / = —AE cosipo, where ipo is the value of i/^ at r = 0, and 



= 2 • u ^ 29 
smh — 

The fact that / displays isolated zeroes as a function of ipo is again proof of the 
presence of chaos when the perturbation is turned on . 



II. 5 - The Separatrix Map 

The oscillatory behavior of the Melnikov Integral proves implicitly that the dynamics 
of our system, restricted to suitable invariant sets, is equivalent to a Bernoulli flow. By 
constructing the so - called "separatrix map''^^^^, we shall be able to show explicitly some 
of those sets, and therefore acquire a much more immediate insight on the evolution of the 
Universe near a resonance. 

To build the separatrix map we shall adopt the picture, emerging from the previous 
subsection, of the dynamics of the conformal field as that of a pendulum, with angle i/j and 
conjugated momentum ^, subject to the perturbation Eq.(27). 

Let us suppose that, at a certain instant r^, the "pendulum" has energy E^. and 
coordinate if^r- If we follow the unperturbed motion over a closed orbit, the energy will 
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keep constant, while the phase (p = t — ip of the perturbation shall increase to (fr+i = 
(Pr + 2T:/ui{Er). Observe that the phase shift is independent of the phase. 

If we consider the effect of the perturbation, over the same lapse, we observe that 
E does not remain constant. For motion close to the separatrix, the change in E is 
approximately the same than at the separatrix itself, and so we find 



Er+i = Er- AE cos (fir (30) 

The energy shift induces changes in the charactheristic frequencies (jj{E) of the un- 
perturbed orbits, and thus the total phase shift is also affected. As a first approximation, 
we may write 

<Pr+l =<Pr + —TE, T (31) 

Eqs. (30) and (31) define the separatrix map. Observe that the phase shift is no 
longer independent of the phase. 

The separatrix map preserves areas. It also has a double sequence of fixed points 
= (£'fc,±7r/2), where k is an integer, and uj{Ek) = l/k. For large k, the fixed points 
approach the separatrix, where^^^] uj{E) ~ ttluq/ \n{32Eo/\Eo — E\), Eq = uj^/Aiv?. Re- 
stricting ourselves to motion below the separatrix, for simplicity, we find Eq — E^ ~ 
32Eq exp{—k7ru>o). 

The behavior of the separatrix map close to a fixed point is determined by its linear 
part. It can be shown that all fixed points are hyperbolic, while the x'j^ points are 
hyperbolic only for large enough k. For our purposes, it is sufficient to concentrate on the 
properties of the map near the x'^ fixed points. 

For large enough k, the eigenvalues of the linearized map around x'j^ are and 1/Afe, 
where 
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Xk^AE—{ 




)\e, 



(32) 



For large fc, Afc ~ {AE /16ujQEo)exp{kTvuJo) diverges exponentially. The map is 
stretching in the direction v'^ — cos(3^d^ + sin (3^dE, where tan/3_|_ ~ AE/{Xk + 1), 
and contracting in the direction = —sinf3-^^p + cos P-ds, where tan /3_ ~ {1/AE){1 — 



We can now see chaos arising. A small parallelogram, with a vertex on xj^ and sides 
along the directions v"*" and v~, is stretched along the former, and contracted along the 
latter. Since v'^ is essentially the (p direction, the angular span of the original parallelogram 
is dilated to many times 27r. After modulo 27r identification, therefore, the image is su- 
perimposed to the original parallelogram, in a way essentially equivalent to Arnold's "Cat 
Map"f^'^^ By a standard procedure, then, it is possible to identify, in our parallelogram, 
an invariant ( null ) set where the separatrix map is Bernoulli t^^. 

Clearly, the procedure works only if the map is strongly dilating. As a ballpark 
estimate, wc can place the "border of chaos" at the value of Ej. for which ~ 1. Thus, 
we are led to estimate the width of the stochastic layer around the departing separatrix 
at 5E ~ 2AE/ujQ. Deep in the resonance region, that is , for the nth resonance, n large, 
Eq.(29) yields AE ~ ujo/inn?' ( recall that ujq ~ n}^^ is large itself ), and so the width of 
the stochastic layer approaches a constant value 5E ~ 2/m^. This value is about four 
times the distance between the resonances themselves, so the stochastic layers merge and 
a stochastic sea is formed: any initial condition in this region falls within the stochastic 
layer of some resonance t^^l 

Almost more important to us, these results show that the loss of information asso- 
ciated to chaotic behavior is observable within the span of a single Universe. Indeed, we 
have identified chaos in the oscillations of the scalar field, parametrized against the angle 
variable associated to the radius of the Universe. Now, over most of the stochastic layer, 
the field performs many oscillations in the lapse between the Big Bang and the Big Crunch. 
This results from the fact that the frequency associated with the expansion and contraction 



(l/Afc)). 
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of the Universe is nearly constant for large m?j, while the frequencies associated to the 
field oscillations scale as ujq ~ n^^^. For example, at the outer rim of the stochastic layer, 
the frequency of field oscillations is ~ (7r/2)(c<;o/ ln2a;o), much larger than the cosmological 
frequency. To reach field frequencies comparable to the cosmic one, we must approach the 
separatrix, up to energies Eq — E 32Eo exp(— 27ra;o). Since the width of the stochastic 
layer does not shrink with n, we see that the fraction of it where the frequency of field 
oscillations is comparable to the rate of cosmic expansion and contraction, is negligible for 
deep resonances. 

To summarize the content of this section, we have shown that for our simple cosmo- 
logical model, the Hamiltonian can be split into an integrable and a non integrable part, 
the integrable one being the main determinant of the dynamics both for small and large 
Universes. Using perturbative methods, we have shown that the perturbation destroys the 
invariant tori of the unperturbed motion. For rn^j ^ 1, j being the adiabatic amplitude 
of field oscillations, a stochastic sea is formed, through the merger of the stochastic layers 
of the individual resonances. Chaos manifest itself through seemingly random phase shifts 
in the field oscillations as the Universe expands and recoUapses. Thus severe limitations 
to our ability to predict the future behavior of the field should arise, even within the span 
of a single Universe. 

Ill - Numerical treatment of the model 

although the analysis in the previous section amounts to rather impressive evidence 
of chaotic behavior in our model, it suffers from the limitations of the perturbative ap- 
proach we have chosen. By solving the model numerically, we shall be able to go beyond 
perturbation theory, and thus find an independent confirmation of the results above. 

In our numerical work, we have used the "physical" variables (a, tt, 0, p), whose equa- 
tions of evolution follow from the Hamiltonian Eq. (5); the simplicity of these equations 
makes this approach more appealing than other, more sophisticated, alternatives. 
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To solve the model, we have used a Runge Kuttat^^l S*'* order routine, implemented 
on an IBM compatible PC 486. As a check on the numerical code, in all runs we have 
surveyed the value of the Hamiltonian, finding that it never exceeded a threshold of 10"^'^. 
For simplicity, we shall set m = 1 throughout. 

Figures (1) and (2) show the evolution of a "typical Universe" in the region of "weak 
chaos" (m^j ~ .5, in the notation of Section II. 1). Figure (1) shows the evolution of the 
scale factor, and Figure (2) that of the field. The nonlinearity is clearly visible, toghether 
with a sharp rise of the amplitude associated to the radius of the Universe in the last 
oscillations. Observe the simultaneous change in the frequency of field oscillations. 

In what follows, we shall subject our model to a series of numerical experiments, 
searching for unambiguous signals of chaotic, rather than complicated, behavior. To this 
end, we shall study its Lyapunov exponents, the projection of an orbit on configuration 
space, and finally the relationship between the values of the field at two consecutive zero 
crossings of the radius of the Universe. 

III.l - Lyapunov exponents. May chaos be quantified? 

It is possible to obtain a quantitative measure of chaos, through the determination 
of Lyapunov exponents. These are defined by considering the deformation of a e-ball 
of initial conditions in phase space, in an ellipsoid of principal axis Aj at time t. If we 
order the axis by size, then the i -th Lyapunov exponent is given by the double limit 
Li — limg^oliiiit^oo((lii Ai)/t). In a Hamiltonian system, the Lyapunov exponents add to 
zero, as a consequence of Liouville theorem; if the system, moreover, is conservative, then 
two of the Lyapunov exponents are zero, reflecting the invariance of the energy shells, and 
the homogeneity of time. Thus, in a conservative system with two degrees of freedom, such 
as ours, only the largest Lyapunov exponent carries meaningful information. 

While a rigorous determination of Lyapunov exponents is usually numerically pro- 
hibitive, it is possible to estimate them, following a proposal by Wolf et ai.^^^l. The idea is 
to construct "local" Lyapunov exponents by taking the average over time of the logarithms 
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of the eigenvalues of the hnearized evolution operator. While these "local" Lyapunov expo- 
nents show a marked transient behavior over short times, they approach the true Lyapunov 
exponent as time goes to infinity. 

Figures (3), (4) and (5) show plots of the largest Lyapunov exponent for three different 
"cosmologies", corresponding to m?j ~ .5, 50 and 5, 000, respectively. In the latter cases, 
we have cut off the early transient, to allow for a larger vertical scale. Beyond the transient 
regime, it can be seen that the greatest Lyapunov exponents are clearly positive, taking 
values in the intervals from .4 to .8, 2 to 3, and 3.5 to 4, respectively. This is consistent 
with the results of Section IL 

III. 2 - Orbits in conGguration space. May chaos be proven? 

In spite of their usefuUness to the determination of chaos, Lyapunov exponents are 
sometimes schewed in studies of relativistic problems, since they are generally gauge depen- 
dent . More reliable information can be obtained from the analysis of Poincare sections. 
A Poincare section is the set of points in which a given orbit of the system intersects, in 
a given direction, a given plane in phase space. If the motion is regular, the Poincare 
sections shall fit in smooth curves, while this will not occur for chaotic motion. 

Similarly, we can study the projection of a trajectory of our system on the (a, 0) 
plane. If a second constant of motion C, in involution with the Hamiltonian H, existed, 
making the system integrablet^^^, then at every point of the trajectory it would be possible 
to solve the equations C — constant and H — constant for n and p, and thereby reduce 
Hamilton's equations to parametric equations for a curve = 4>{a). For example, if m = 0, 
then we may take C = {p^ + 0^)/2. In these circumstances, the projection of the orbit on 
the (a, 0) plane fits into a smooth curve. Faillure to do so, on the other hand, evidences 
the non existence of a second constant of motion ( disregarding pathological cases ). 

Figures (6), (7) and (8) display plots of vs a, again for m^j equal to .5, 50 and 
5,000, respectively. It can be seen that, while in the first case it seems possible to fit the 
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resulting graph into a smooth Lissajous curve, this is no longer so in the latter cases. This 
again confirms the results from Section II. 

III. 3 - Field values at consecutive cosmic singularities. May chaos be observed? 

In the above numerical tests, we have allowed the radius of the Universe to cross zero 
several times. However, in the "real" world our observations cannot be projected beyond 
the cosmological singularity. Thus it is crucial to determine whether the loss of accuracy 
in the prediction of the future behavior of orbits, associated to chaos, is strong enough to 
lead to observable results within a single Universe. 

In order to do that, we analyzed what happens with Universes starting from the same 
geometrical condition ( that is, same a and tt ) but different values of </> ( and p given by 
the Hamiltonian constraint ). We begin our computation a little while after the Big Bang, 
ending the numerical calculation when a becomes negative for the first time ( that is, on 
the Big Crunch ). Repeating this procedure several times, a plot of the final value 0/ of 
the field against the initial one 0j can be obtained. For regular behavior, this plot should 
allow the prediction of (f)f given 0j with the proper accuracy. 

Figures (9), (10) and (11) show the results of numerical experiments such as described, 
once more for m^j ~ .5, 50 and 5,000, respectively. In the plot Fig. (9), all simulations 
started at a = .001 and tt = —1, while 0j ranged from .001 to 1; to obtain Fig (10), we set 
the initial values of a and tt at .001 and —10, respectively, allowing to range from .01 
to 10; finally. Fig. (11) corresponds to the values a, = .1, Tij = —100, .1 < 0i < 100. 

The outcome of the simulations indicates that, while 0/ is easily predicted from 0j 
in the first case, prediction becomes impossible in the latter ones. Indeed, the correlation 
coefficient between (pf and (pi, which gives .933 for the data in Fig. (9), falls to .049 and 
.103 for Figs. (10) and (11), respectively. Moreover, regular "islands" which can still be 
seen in Fig (10) ( e. g., for 0j between .3 and .45 ) have quite dissapeared from Fig (11). 

In conclusion, numerical simulation of the model strongly confirms the results of the 
analysis in Section II. Moreover, chaos places severe limitations in our ability to predict 
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the future behavior of the Universe, even before the Cosmic Singularity is reached for the 
first time. 

IV - Conclusions 

The analytical and numerical analysis above shows that the dynamics of a spatially 
closed Friedmann - Robertson - Walker Universe coupled to a conformal but massive scalar 
field displays chaotic behavior. Indeed, in appropiate invariant sets, the dynamics of this 
cosmological model is equivalent to a Bernoulli shift, and thus essentially indistinguishable 
from a purely random process^^*^!. The analysis above tells us how to identify those sets. 

The random element we have alluded to appears in the solutions of the model through 
the nonlinear oscillations of the field as the Universe evolves from the initial cosmic sin- 
gularity towards recoUapse. As a result, our ability to predict the future behavior of the 
field is limited. This phenomenon demonstrably occurs within a single cosmic episode, and 
thus it is, in principle, observable. 

The interest of this result probably lies less in its direct impact on present theories 
of the Universe, than in the avenues it opens for further research. As a description of the 
Universe, the model we have discussed here is obviously oversimplified. More seriously, 
the chaotic behavior we have analyzed requires high values of the field energy density, a 
regime where quantum effects should not be disregarded. 

In our view, the interest of this model lies in its sheer simplicity. In this sense, we 
believe this model to be an almost ideal ground where to investigate the general issues 
associated with chaos and cosmology. For one thing, this model is so simple that its 
analysis at the quantum and semiclassical levels is not essentially harder than at the 
classical one. This affords an excellent opportunity to discuss the semiclassical limit of 
Quantum Cosmology in a highly non trivial framework, an issue which we will discuss in 
a separate communication. 
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On the other hand, the model we have presented is far from displaying the full richness 
of behavior than chaos may bring to cosmology. We continue our research on this most 
varied and rewarding field. 
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FIGURE CAPTIONS 



Figure 1 Evolution of ttie scale factor for a Robertson- Walker cosmology in a weak chaotic 
situation (^j = 1 and tTj = —1). After several oscilations there is an unusual growth of the 
scale factor, this corresponds to an "inflationary-like" stage. 

Figure 2 Evolution of the field (p corresponding to the inital conditions of Fig.l. It can 
be appreciated the oscillatory behavior with a sudden change, in frecuency and amplitude, 
during the "inflationary " stages. 

Figure 3 Maximal Lyapimov exponent {(pi = 1, iVi — —1 and discarding a transient stage 
near rj = 0). Its value is clearly greater than zero. 

Figure 4 Maximal Lyapunov exponent {(pi = 10, TTj = —10 and discarding a transient 
stage near r; = 0). 

Figure 5 Maximal Lyapunov exponent {(pi = 1, TTj = —100 and discarding a transient 
stage near 77 = 0). It can be appreciated that its value increase with the increase of initial 
condition, as predicted in section II. 

Figure 6 Projection of the phase space trajectory onto the (p , a plane. In this situation 
{(pi = 1 and TTi = —1), it behaves as a smoth Lissajous curve. 

Figure 7 Same as Fig. 6, {(pi = 1 and tt^ = —100). The trajectory is very irregular, as 
indicated by its high Lyapunov exponent. 

Figure 8 Same as Fig. 6. This corresponds to the Lyapunov shown in Fig.4 {(pi — 10 and 
TTi = —10). See the increase of amplitude in the scale factor in the last two cases. 

Figure 9 Field at the Big Crunch against field at the Big Bang. The initial data are: 
tti = 0.001 (for computational conveniences), tTj = — 1 and (pi running from zero to one. 
There is a simple relation between (pf and (pi. 

Figure 10 Same as Fig. 9, with = 0.001, tt^ = —10 and (pi running from zero to ten. 
"Islands" of stability may still be seen, for example between (1.5 and 3), (3 and 4.5) and 
(5 and 6). 
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Figure 11 Same as Fig. 9, with = 0.1, tTj = —100 and 0j running from zero to one 
hundred. In this situation there is no clear correlation between initial and final data. 
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